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Abstract. In this article we present a new class of high order accurate Arbitrary- 

Eulerian-Lagrangian (ALE) one-step WENO finite volume schemes for solving non- 
linear h5^erbolic systems of conservation laws on moving two dimensional unstruc- 
tured triangular meshes. A WENO reconstruction algorithm is used to achieve high 
order accuracy in space and a high order one-step time discretization is achieved by 
using the local space-time Galerkin predictor proposed in [25]. For that purpose, a 
new element-local weak formulation of the governing PDE is adopted on moving 
space-time elements. The space-time basis and test functions are obtained consid- 
ering Lagrange interpolation pol3momials passing through a predefined set of nodes. 
Moreover, a pol5momial mapping defined by the same local space-time basis func- 
tions as the weak solution of the PDE is used to map the moving physical space-time 
element onto a space-time reference element. To maintain algorithmic simplicity, the 
final ALE one-step finite volume scheme uses moving triangular meshes with straight 
edges. This is possible in the ALE framework, which allows a local mesh velocity that 
is different from the local fluid velocity. We present numerical convergence rates for 
the schemes presented in this paper up to sixth order of accuracy in space and time and 
show some classical numerical test problems for the two-dimensional Euler equations 
of compressible gas dynamics. 
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1 Introduction 

In this paper we present a new family of high order accurate Lagrangian-type one-step 
finite volume schemes for solving nonlinear hyperbolic balance laws, with non stiff al- 
gebraic source term. The main advantage of working in a Lagrangian framework is that 
material interfaces can be identified and located precisely, since the computational mesh 



* Corresponding author. Email addresses: walter.boscheri@imitn.it (W. Boscheri), 
michael . dimbserOiinitn . it (M. Dumbser) 



http://www.global-sci.com/ 



Global Science Preprint 



2 



is moving with the local fluid velocity, hence obtaining a more accurate resolution of ma- 
terial interfaces. For this reason a lot of research has been carried out in the last decades 
in order to develop Lagrangian methods, whose algorithms can start either directly from 
the conservative quantities such as mass, momentum and total energy [57,68], or from the 
nonconservative form of the governing equations, as proposed in [6,9,79]. Furthermore 
we can split the existing Lagrangian schemes into two main classes, depending on the lo- 
cation of the physical variables on the mesh: in one case the velocity is defined at the cell 
interfaces and the other variables at the cell barycenter, hence adopting a staggered mesh 
approach, while in the other case all variables are defined at the cell barycenter, therefore 
using a cell-centered approach. 

The equations of Lagrangian gas dynamics have been considered in [60], where sev- 
eral different Godunov-type finite volume schemes have been presented and where a 
new Roe linearization has been introduced in order to define proper estimates of the 
maximum signal speeds in HLL-type Riemann solvers in Lagrangian coordinates. A 
cell-centered Godunov scheme has been proposed by Carre et al. [10] for Lagrangian gas 
dynamics on general multi-dimensional unstructured meshes. In this case the finite vol- 
ume scheme is node based and compatible with the mesh displacement. In [21] Despres 
and Mazeran introduce a new formulation of the multidimensional Euler equations in La- 
grangian coordinates as a system of conservation laws associated with constraints. Fur- 
thermore they propose a way to evolve in a coupled manner both the physical and the ge- 
ometrical part of the system [22], writing the two-dimensional equations of gas dynamics 
in Lagrangian coordinates together with the evolution of the geometry as a weakly hy- 
perbolic system of conservation laws. This allows the authors to design a finite volume 
scheme for the discretization of Lagrangian gas dynamics on moving meshes, based on 
the symmetrization of the formulation of the physical part. In a recent work Despres et 
al. [18] propose a new method designed for cell-centered Lagrangian schemes, which is 
translation invariant and suitable for curved meshes. General polygonal grids are also 
considered by Maire et al. [54-56], who develop a general formalism to derive first and 
second order cell-centered Lagrangian schemes in multiple space dimensions. By the use 
of a node-centered solver [56], the authors obtain the time derivatives of the fluxes. The 
solver may be considered as a multi-dimensional extension of the Generalized Riemann 
problem methodology introduced by Ben-Artzi and Falcovitz [5], Le Floch et al. [8, 41] 
and Titarev and Toro [70, 71, 73]. So far, all the above-mentioned schemes are at most 
second order accurate in space and time. 

In order to achieve higher accuracy, Cheng and Shu were the first who introduced 
a high order essentially non-oscillatory (ENO) reconstruction in Lagrangian schemes 
[14,53]. They developed a class of cell centered Lagrangian finite voliune schemes for 
solving the Euler equations, using both Runge-Kutta and Lax-Wendroff-type time step- 
ping to achieve also higher order in time. Furthermore a formalism for symmetry pre- 
serving Lagrangian schemes has been proposed by Cheng and Shu, see [15, 16]. The 
higher order schemes presented in [14, 53] were the first better than second order non- 
oscillatory Lagrangian-type finite voliuxie schemes ever proposed. Higher order finite 
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element methods on unstructured meshes have been investigated in [62, 67]. In a very 
recent paper Dumbser et al. [30] propose a new class of high order accurate Lagrangian- 
type one-step WENO finite volume schemes for the solution of stiff hyperbolic balance 
laws. They consider the one-dimensional case and develop a Lagrangian algorithm up 
to eighth order of accuracy, based on high order WENO reconstruction in space and the 
local space-time discontinuous Galerkin predictor method proposed in [26] to obtain a 
high order one-step scheme in time. 

Other schemes can be also included and mentioned as Lagrangian algorithms, e.g. 
meshless particle schemes that adopt a fully Lagrangian approach, such as the smooth 
particle hydrodjmamics (SPH) method [36-39, 58], which can be used to simulate fluid 
motion in complex deforming domains. Also within the SPH approach, which is a fully 
Lagrangian method, the mesh moves with the local fluid velocity, whereas in Arbitrary 
Lagrangian Eulerian (ALE) schemes, see e.g. [13, 23, 34, 35, 44, 64, 68], the mesh moves 
with an arbitrary mesh velocity that does not necessarily coincide with the real fluid 
velocity. Furthermore one can find Semi-Lagrangian schemes, which are mainly used for 
solving transport equations [42,66]. Here, the numerical solution at the new time level 
is computed from the known solution at the present time by following backward in time 
the Lagrangian trajectories of the fluid to the end-point of the trajectory. Since the end- 
point does not coincide with a grid point, an interpolation formula is required in order 
to evaluate the imknown solution, see e.g. [7,11,12,20,46,52,65]. For the sake of clarity 
we specify that in Semi-Lagrangian algorithms the mesh is fixed as in a classical Eulerian 
approach. An alternative to Lagrangian methods for the accurate resolution of material 
interfaces has been developed in the Eulerian framework on fixed meshes under the form 
of the ghost-fluid method [32,33,40], together with a level-set approach [59,63], where 
the level-set function represents the signed distance from the material interface and its 
zeros locate the interface position. 

In this paper we introduce a new better than second order accurate two-dimensional 
Lagrangian-type one-step WENO finite volume scheme on unstructured triangular meshes: 
high order of accuracy in space is obtained using a WENO reconstruction [2, 3, 27, 28, 
45, 47, 72, 76, 80], although other higher order spatial reconstruction schemes could be 
adopted as well, see [1, 17]. High order accuracy in time is achieved with a local space- 
time Galerkin predictor, as introduced in [26,30,31,43]. The method proposed in this arti- 
cle is presented as an Arbitrary-Lagrangian-Eulerian scheme in order to allow arbitrary 
grid motion. This allows us to use curved space-time elements in the local predictor stage 
but triangles delimited by straight line segments in the resulting one-step finite volume 
scheme (corrector stage). This choice has been made to maintain algorithmic simplicity. 

The outline of this article is as follows: in Section 2 we present the details of the 
proposed niunerical scheme, while in Section 3 we show numerical convergence rates up 
to sixth order of accuracy in space and time for a smooth problem as well as numerical 
results for several different test cases governed by the compressible Euler equations. The 
paper closes with some concluding remarks and an outlook to possible future extensions 
of the method in Section 4. 
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2 Numerical Method 

In this article we consider general nonlinear systems of hyperbolic balance laws of the 
form 

^+V-F(Q) = S(Q), (x,y)eO(0cR2, ieR+ Qe^QC]R^ (2.1) 

where Q = {(]\,q2,-;qv) is the vector of conserved variables defined in the space of the 
admissible states Oq C W , F(Q) = (f (Q),g(Q)) is the nonlinear flux tensor and S(Q) 
represents a nonlinear but non-stiff algebraic source term. The spatial position vector is 
denoted by x= {x,y), t is the time and Q(t) is the time-dependent computational domain. 

The two-dimensional time-dependent computational domain Cl{t) is discretized at 
the current time t" by a set of triangular elements T". The union of all elements is called 
the current triangulation Tq of the domain 0(i") =0" and can be expressed as 

T^ = [jTf, (2.2) 

i=l 

where Ne is the total number of elements contained in the domain. 

In a Lagrangian framework we are dealing with a moving mesh, hence the elements 
deform while the solution is evolving in time. It is therefore convenient to adopt a local 
reference coordinate system £,—t}, where the reference element Tg is defined. The spatial 
mapping of the triangular elements T" at the current time P from reference coordinates 
f — ?/ to physical coordinates x— y is given by the relation 

y = yi,i+{yv-yi,i)^+{ys,-yi,i)v^ (2.3) 

where • = denotes the vector of spatial coordinates of the A;-th vertex of the 

triangle T" in physical coordinates at the current time t". The vector ^={^,tj) is the vector 
of spatial coordinates in the reference system, while x=(x,y) is the spatial coordinate vec- 
tor in the physical system. The spatial reference element Tg is the unit triangle composed 
of the nodes ^ = (l,c,\,r]e,i) = (0,0), 2 = i^e,2,ne,2) = (1,0) and 3 = {l,e,2,r]e,2) = (0,1). 

As usual for finite volume schemes, data are represented by spatial cell averages, 
which are defined at time P as 

Q" = ^ f Qix,yr)dV, (2.4) 

where \T"\ denotes the volume of element T" at the current time t". To achieve higher 
order in space, piecewise high order polynomials w;,(x,y,i") must be reconstructed from 
the given cell averages (2.4) using a pol5niomial WENO reconstruction procedure illus- 
trated briefly in the next section. 
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2.1 Polynomial WENO Reconstruction on Unstructured Meshes 

In this paper we use the WENO reconstruction algorithm in the polynomial formulation 
presented in [26-28], instead of adopting the original pointwise WENO scheme [45,47,80]. 
Here we will give a very brief summary of the algorithm, since it is completely described 
in all detail in the above-mentioned references. 

The reconstruction is performed by using the reference system {^,t]) according to the 
mapping (2.3), as also explained in detail in [27]. The reconstruction polynomial of degree 
M is obtained by considering a reconstruction stencil 

S! = [jT:(jy (2.5) 

;=i 

The stencil contains a total number of elements ng that is greater than the smallest number 
M. = (M+l)(M+2)/2 needed to reach the formal order of accuracy M+1, as shown 
in [4,49,61]. Typically we take ne=2^A in two space dimensions. Here 1 <;' < is a local 
index counting the elements belonging to the stencil, while m{)) maps the local index to 
the global element numbers used in the triangulation (2.2). s denotes the number and 
the position of the stencil with respect to the central element T". According to [27,28] 
we always will use one central reconstruction stencil given by s = 0, three primary sector 
stencils and three reverse sector stencil, as suggested by Kaser and Iske [49], so that we 
globally deal with seven stencils per element. 

We use some spatial basis functions tpi{^,t]) in order to write the reconstruction poly- 
nomial for each candidate stencil s for triangle T": 

M 

<{^>y>n = EM^'V)K: ■.=m,fj)ikl'!> (2-6) 

i=\ 

with the mapping to the reference coordinate system given by the transformation (2.3). 
In the rest of the paper we will use classical tensor index notation with the Einstein sum- 
mation convention, which implies summation over two equal indices. The number of the 
unknown polynomial coefficients (degrees of freedom) is M., already defined previously. 
As basis functions tpi{^,r]) we use the orthogonal Dubiner-tj^pe basis on the reference 
triangle Tg, given in detail in [19,24,48]. 

Integral conservation is required for the reconstruction on each element T"eS-, hence 

^|V'/(e,'/)wr;W=Q«, ^TfGSf (2.7) 

where |T"| represents the voliune of element T" at time P. Since the number of sten- 
cil elements is larger than the one of the unknown polynomial coefficients (n^ > M), the 
above system given by eqn. (2.7) is an overdetermined linear algebraic system that is 
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readily solved for the unknown coefficients w'^f using a constrained least-squares tech- 
nique, see [27]. The multi-dimensional integrals are evaluated using Gaussian quadra- 
ture formulae of suitable order, see the book of Stroud for details [69]. Since the shape 
of the triangles changes in time, the small linear systems given by (2.7) must be solved 
at the beginning of each time step, while the choice of the stencils S- remains fixed for 
all times. It is therefore very useful to devise a one-step time-discretization, such as the 
one detailed in the next section, which requires only one reconstruction per time step, in 
contrast to higher order Runge-Kutta methods, which would require the evaluation of 
the above reconstruction operator in each substage of the Runge-Kutta method. 

In order to obtain essentially non-oscillatory properties the final WENO reconstruc- 
tion polynomial is computed from the reconstruction polynomials obtained on each in- 
dividual stencil S? in the usual way. For this purpose we adopt the classical definitions 
of the oscillation indicators os reported in [47] and the oscillation indicator matrbc Z/^ 
proposed in [27,28], which read 



The nonlinear weights tOg are defined by 



As a;s ,^ „ 

= 7 — —Tf, cOs = =^-^, (2.9) 

where we use e = 10"^^, r = 8, As = 1 for the one-sided stencils and Aq = 10^ for the central 
stencil, according to [26,28]. The final nonlinear WENO reconstruction polynomial and 
its coefficients are then given by 



M 

vfh{x,y,n = Y,%{^,n)^li> with wr, = E^swff. (2.10) 

1=1 



2.2 Local Space-Time Continuous-Galerkin Predictor on Moving Meshes 

In order to obtain high order accuracy in time, the reconstructed polynomials W;, ob- 
tained at the current time P are now evolved locally within each element T, (f) during one 
time step t" <t< t"^^. The result of this local evolution step is an element-local piecewise 
polynomial space-time predictor qh{x,y,t). For this purpose we use an element-local 
weak formulation of the PDF in both space and time, based on a Lagrangian version of 
the local space-time continuous Galerkin method introduced for the Fulerian framework 
in [25]. 

Let X = {x,y,t) denote the physical coordinate vector and f = (^,//,t) the reference co- 
ordinate vector, in both of which also the time is included, while x = {x,y) and ^ = (^,?/) 
are the pure spatial coordinate vectors in physical and reference coordinates, respectively. 
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Let furthermore 0/ = = 9i{^,f],T) be a space-time basis function defined by the La- 
grange interpolation polynomials passing through the space-time nodes = (^,„,//,„,Tm). 
The nodes are specified according to [25]. Being the Lagrange interpolation polynomials, 
the resulting nodal basis functions 9i therefore satisfy the interpolation property 

e,{U=^im. (2.11) 

where 3i„ is the usual Kronecker symbol. According to [25] the local space-time solution 
in element Ti{t), denoted by q^, as well as the fluxes (f;„g;,) and the source term S;, are 
approximated as 

h = hi^,rj,T) = 01 {^,r],T)\i, g/, = g/, {^,r],T) = Oi (^,7/,T)g/,,. (2.12) 

The mapping between the physical space-time coordinate vector x and the reference 
space-time coordinate vector f is achieved using an isoparametric approach, i.e. the map- 
ping is represented by the same basis functions 9i as the approximate niunerical solution 
itself, given in Eqn. (2.12) above. Hence, one has 

x(^,?/,t)=0/(^,?/,t)%, y{^,r],r) = 6i{^,r],r)yii, i(^,;/,T) =0/(^,?/,t)T/, (2.13) 

where the degrees of freedom x/ /= (£;„y/ ;) denote the (in parts unknown) vector of phys- 
ical coordinates in space of the moving space-time control volume and the ti denote the 
known degrees of freedom of the physical time at each space-time node = {xii,yii,ti). 
The mapping in time simply reads 

t-t" 

t = tn+rAt, '^=-^' =^ ti = tn+riM, (2.14) 

where is the current time and M is the time step. Hence, t^ = trj=0 and tx = At. 

The spatial mapping given in (2.13) and the temporal mapping (2.14) allow us to 
transform the physical space-time element to the unit reference space-time element Tg x 
[0,1]. The Jacobian of the spatial and temporal transformation is given by 

(2.15) 




and its inverse reads 




(2.16) 



The local reference system and the inverse of the associated Jacobian matrix (2.16) are 
used to rewrite the governing PDE (2.1) as 
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which simplifies to 

r;^n ?in ?i( 7)( 7)0 1 

:AfS(Q), (2.18) 



^ '-At 



9Q^ dQ af^ 3f 9g^ 9g 



since = Ty = and Tf = according to the definition (2.14). 

To make the notation easier we now introduce the following two integral operators 

1 

[f'gV = I mri'T)g{^,f]'r)d^df], {f,g) f{^,f],T)g{^,fj,T)d^dridT, (2.19) 

Te % 

which denote the scalar products of two functions / and g over the spatial reference 
element Tg at time r and over the space-time reference element Tg x [0,1], respectively. 

Inserting the definitions of (2.12) into the weak formulation of the PDE (2.18), we mul- 
tiply Eqn. (2.18) with the same space-time basis functions 6]^{^,}],t) and then integrate it 
over the space-time reference element Tg x [0,1], hence obtaining: 

+ At (^^% ^^y^ + (%^'?y)) hi = M) Si,i- 

(2.20) 

The above expression can be written in a more compact matrix form, which reads 

Krqi,i+At (Ktqi,i+KJi^i+Ky%i^^, = AtM%, (2.21) 
with the following matrix definitions: 

Kr = /%|^V M = {ekA), (2.22) 



dr 

and 
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Note that the Lagrangian nature of the scheme, i.e. the moving space-time control vol- 
ume, leads to the term Ktqu, which is not present in the Eulerian case introduced in [25]. 
In (2.21), only the matrix Kt and the space-time mass matrix M given in (2.22) can be 
pre-computed and stored. The other integrals continuously change in time, since they 
depend on the inverse of the Jacobian matrix, which varies according to the Lagrangian 
mesh motion. Therefore it may be more convenient to assemble a unified term P as 

hence 

Qr = MF. (2.25) 

P is approximated again inside element T,(i) by adopting a nodal approach as done in 
(2.12), i.e. 

Fh = n{^,V>^)=H^>V>^)hu (2.26) 

withP,, = P(x,,). 

Let us denote with qj*- the degrees of freedom that are known from the initial con- 
dition w;, so that (\h{x,y,t") =vf}i{x,y,t"), whereas the unknown degrees of freedom for 
T > are denoted by qj^, so that the total vector of degrees of freedom is written as 
qi i = (q°,,q^,). One can move onto the right-hand side of (2.21) or (2.25) all the known 
degrees of freedom q^- given by the initial condition w;, by setting the corresponding 
degrees of freedom to the known values, like a standard Dirichlet boundary condition 
in the continuous finite element framework, see [25] for more details. The final expres- 
sion for the iterative scheme solving the nonlinear algebraic equation system of the two- 
dimensional Lagrangian continuous Galerkin predictor method simply reads 

Krq[f = AiMP[,., (2.27) 

where the superscript r denotes the iteration number here. For an efficient initial guess 
(r = 0) based on a second order MUSCL-type scheme, see [43]. 

Since the mesh is moving we also have to consider the evolution of the vertex coordi- 
nates of the local space-time element, whose motion is described by the ODE system 

d\ 

^=V(x,y,f), (2.28) 

where V=V(x,i/,f) = {U,V) is the local mesh velocity. Since we are developing an arbitrary 
Lagrangian-Eulerian scheme (ALE), we want the mesh velocity to be independent from 
the fluid velocity. In this framework we can deal in the same way with both Eulerian and 
Lagrangian schemes: if V = the scheme reduces indeed to a pure Eulerian approach, 
while if V coincides with the local fluid velocity v we obtain a Lagrangian-type method. 
The velocity inside element Tj(f) can be also expressed as 



(2.29) 
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withV,,=V(x/,;). 

As introduced in [30], we can solve the ODE (2.28) for the unknown coordinate vector 
X/ using again the local space-time CG method: 



hence resulting in the following iteration scheme for the element-local space-time pre- 
dictor for the nodal coordinates: 



The nodal degrees of freedom x/ at relative time t=0, i.e. the initial condition of the ODE 
system, are given by the spatial mapping (2.3), since the physical triangle T" at time t" is 
known. 

Eqn. (2.30) is iterated together with the weak formulation for the solution (^.17). The 
algorithm for the two-dimensional local space-time CG predictor for moving meshes can 
be summarized by the following steps: 

• compute the local mesh velocity (2.29), usually by choosing the fluid velocity, hence 

Vu=v(xm); 

• with (2.31) update the geometry locally within the predictor stage, i.e. the element- 
local space-time coordinates; 

• compute the Jacobian matrix and its inverse by using (2.15)-(2.16); 

• compute the term P according to (2.25); 

• evolve the solution locally with (2.27). 

The iterative procedure described above stops when the residuals of (2.27) and (2.31) 
are less than a prescribed tolerance tol (typically tol^lQ~^^). 

Once we have carried out the above procedure for all the elements of the computa- 
tional domain, we end up with an element-local predictor for the numerical solution q;,, for 
the fluxes F;, = (f;„g;,), for the source term S;, and also for the mesh velocity V;,. 

Next, we have to update the mesh globally. Let us denote with Vk the neighborhood 
of vertex number k, i.e. all those elements that have in common the node nimiber k. The 
number of elements in the neighborhood is denoted with N^. Since the velocity of 
each vertex is defined by the local predictor within each element, one has to deal with 
several, in general different, velocities for the same node, since all elements belonging to 
Vfc will in general give a different velocity contribution, according to their element-local 
predictor. Since we do not admit the geometry to be discontinuous, we decide to fix a 




(2.30) 



(2.31) 
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fl 

unique node velocity V^^ to move the node. The final velocity is chosen to be the average 
velocity considering all the contributions Vjjy of the vertex neighborhood as 

Vfc = ^ %i' with Vlj =(^l9i iUmikyVeMk)'^)'^^'^ ^H- (2-32) 

The ^e,m{k) ^i^d t}e,m{k) ^^e the vertex coordinates of the reference triangle Tg corresponding 
to vertex number k, hence m = m{k) with 1 < m < 3 is a mapping from the global node 
number k to the element-local vertex number. Since each node now has its own imique 
velocity, the vertex coordinates can be moved according to 

X^+i=X^+AiV^, (2.33) 

and we can update all the other geometric quantities needed for the computation, e.g. 
normal vectors, volumes, side lengths, barycenter position, etc. 

2.3 Finite Volume Scheme 

The conservation law (2.1) can be easily rewritten in a more compact space-time diver- 
gence form, which reads 

V-Q = S(Q), (2.34) 

with 



Integration over the space-time control volume C" = Ti{t) x yields 

f«+l jn+l 

J J V-Qdxdt= J J Sdxdt. (2.36) 

f" T,(f) t" T;(t) 

Applying Gauss' theorem allows us to write the left space-time volume integral above 
as sum of the flux integrals computed over the space-time surface dC", hence 



Q-hdS= J J Sdxdt, (2.37) 

dCf t" Ti(t) 

where h={nx,ny,nt) is the outward pointing space-time unit normal vector on the space- 
time surface dCf. 

The space-time surface dCf above involves overall five space-time sub-surfaces, as 
depicted in Figure 1: 

dcf = I IJ ac;;. I u Tf u Tr'+\ (2.38) 
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where Mi denotes the so-called Neumann neighborhood of triangle Ti{t), i.e. the set of 
directly adjacent triangles Tj{t) that share a common edge dTij{t) with triangle Ti{t). 
The common space-time edge dC"- during the time interval [t";t"^^] is denoted above 

by8q=aTy(0x[r;r+i]. 

The upper space-time sub-surface T"^^ and the lower space-time sub-surface T" are 
parametrized byO<f<lAO<?/<l — ^ and the mapping (2.3). They are orthogonal to 
the time coordinate, hence for these faces the space-time unit normal vectors simply read 
ft = (0,0,1) for Tf+^ and ft = (0,0,-1) for Tf, respectively. 

The lateral space-time sub-faces dCfj are defined using a simple bilinear parametriza- 

tion, since the old vertex coordinates X"^ are given and the new ones X"^"^^ are known from 
(2.33). 

dCfj = x{x,T) = '£^k{Xrr)Xl„ 0<X<h 0<T<1, (2.39) 

fc=i 

where (;y,t) represents a side-aligned local reference system according to Figure 1. The 
Xfj J, are the physical space-time coordinate vectors for the four vertices that define the 
lateral space-time sub-surface dC"j. If X"j j and X"y2 denote the two spatial nodes at time 
P that define the common spatial edge dTij{t"), then the four vectors X| ^ are given by 

^>n ( -yn i.n\ -yn / -yn in\ -yn / -yn+l i.n+l\ v-n / -yn+l i.n+l\ 

J' ^ii,2—\^ii,2>^ J> ^ii,3— \^ij,2 >^ )' ^ijA~ )■ 

(2.40) 

The ^kiX''^) ^re a set of bilinear basis functions, which are defined as 

^i{x,r) = {l-x){l-r), 

hi.x,'^) = x^, 

^^ix,r) = il-x)r. (2.41) 

From (2.40) and (2.41) it follows that the temporal mapping is again simply t = t"+rAt, 
hence t^ = and tr = M. 

The determinant of the coordinate transformation and the resulting space-time imit 
normal vector ft,y of the sub-surface dC"- can be computed as follows: 



laqi 



Therefore, the ALE-type one-step finite volume scheme takes the following form: 

1 1 t"+i 
|^„+i|Q„+i^|^„|Qn_ ^ I l\dCf.\Qij.hijdxdr+ J Js{qn)dxdt, (2.43) 



TjeAfi'Q t" T,(t) 




where Q;y • n,y is an Arbitrary-Lagrangian-Eulerian numerical flux function to resolve the 
discontinuity of the predictor solution qjj at the space-time sub-face dC"j. The surface in- 
tegrals appearing in (2.43) are approximated using multidimensional Gaussian quadra- 
ture rules, see [69] for details. At the interface 9C-j let us denote the local space-time 

predictor solution inside element T;(t) by q^^^^ and the element-local predictor solution of 
the neighbor element Tj{t) by q^, then a simple ALE-type Rusanov flux [30] is given by 

Q,7-fiy = ^(Q(q,|)+Q(q,7))-fiy-^s^ax(q,|-q,-), (2.44) 

where Smax is the maximum eigenvalue of the ALE Jacobian matrix in spatial normal 
direction, 

A^:(Q) = 9(F-n)/3Q-(V-n)I, with n = -^^^^f^. (2.45) 

Here, V ■ n is the local normal mesh velocity and I is the vxv identity matrix. It can be 
easily verified that 

V= — f^^Y hu=( -IIm I, hence V-n = -^^=. (2.46) 

A more sophisticated Osher-type ALE flux has been introduced in the Eulerian case 
in [29] and has been employed in the Lagrangian framework in [30]. It reads 

Q,yn,-, = ^(Q(q+) + Q(q, ))-n,,-^ [ J | aX(Y(s)) | j (q,+ -q,), (2.47) 
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with the straight-line segment path connecting the left and right state as follows, 

Y(s) = q,-+s(q+-q,-), 0<s<l. (2.48) 

In (2.47) above, the usual definition of the matrix absolute value operator applies, i.e. 

|A|=R|A|R-i, |A|=diag(|Ai|,|A2|,...,|A^|), (2.49) 

with the right eigenvector matrix R and its inverse R^^. According to [29] the path in- 
tegral appearing in (2.47) is approximated using Gaussian quadrature rules of sufficient 
accuracy. 

3 Test Problems 

In order to validate the two-dimensional ALE finite volume scheme presented in the 
previous section, we present some computational test problems in the following. We 
consider the Euler equations of compressible gas dynamics, which read: 

Qt+fx+gy = Six,y,t) (3.1) 

with 



/ p \ 




1 pu \ 




/ pv \ 


pu 

pv 


, f= 


pu'^ + p 

puv 


' g = 


puv 

pv^ + p 


\ PE ) 




\ u{pE + p) ) 




\ ^{pE + p) 1 



(3.2) 



Here, p denotes the fluid density , v= {u,v) is the fluid's velocity vector and pE represents 
the total energy density. The vector of source terms is denoted by S and the fluid pressure 
by p, given in terms of the conserved quantities by the equation of state of an ideal gas as 

p = {^-l)(^pE-^-p{u^W)y (3.3) 

with the ratio of specific heats 7. 

For each of the following test cases we choose the local mesh velocity as the local fluid 
velocity according to Eqn. (2.32), i.e. 

V = v. (3.4) 
3.1 Numerical Convergence Results 

The numerical convergence studies for the two-dimensional Lagrangian finite volume 
scheme are carried out considering a smooth convected isentropic vortex, see [45]. The 
initial condition is given in terms of primitive variables and it consists in a linear super- 
position of a homogeneous background field and some perturbations 5: 



{p,u,v,p) = {l+6p,l+5u,l+5v,l+6p). 



(3.5) 
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Table 1: Numerical convergence results for the compressible Euler equations using the first up to sixth order 
version of the two-dimensional Lagrangian one-step WENO finite volume schemes presented in this article. 
The error norms refer to the variable p (density) at time t = 1.0. 



h{n(tf)) 


^L2 


0{L2) 


h{n,tf) 




0{L2) 


h{n,tf) 




0{L2) 




Ol 






02 






03 




3.73E-01 


9.525E-02 




3.43E-01 


1.716E-02 




3.28E-01 


1.614E-02 




2.63E-01 


6.907E-02 


0.9 


2.49E-01 


1.109E-02 


1.4 


2.51E-01 


6.943E-03 


3.0 


2.14E-01 


5.700E-02 


0.9 


1.69E-01 


5.766E-03 


1.7 


1.68E-01 


2.290E-03 


2.7 


1.74E-01 


4.752E-02 


0.9 


1.28E-01 


3.027E-03 


2.3 


1.28E-01 


9.274E-04 


3.3 




04 






05 






06 




3.29E-01 


4.717E-03 




3.29E-01 


4.9463E-03 




3.29E-01 


2.051E-03 




2.51E-01 


1.822E-03 


3.5 


2.51E-01 


1.4648E-03 


4.5 


2.51E-01 


5.803E-04 


4.7 


1.67E-01 


4.379E-04 


3.5 


1.67E-01 


2.5937E-04 


4.3 


1.67E-01 


8.317E-05 


4.8 


1.28E-01 


1.313E-04 


4.4 


1.28E-01 


6.9664E-05 


4.9 


1.31E-01 


1.994E-05 


5.9 



We set the vortex radius = {x — 5)^-{-{y — 5)^, the vortex strength e = 5 and the ratio of 
specific heats 7 = 1.4. The perturbation of entropy S = ^ is assumed to be zero, while the 
perturbations of temperature T and velocity v are given by 

«=0, n^- '-l-'f e^-. (3.6) 

From (3.6) it follows that the perturbations for density and pressure are given by 

^p=(l+^T)^-l, ^p = (l+^T)^-l. (3.7) 

The initial computational domain Q(0) = [0;10] x [0;10] is square-shaped and it is 
surrounded by four periodic boundary conditions. The exact solution Qg is the time- 
shifted initial condition given by Qe(x,i) = Q(x— Vci,0), with the convective mean veloc- 
ity v^ = (1,1)- We use the Osher-type numerical flux (2.47) to run this test case until a 
final time of = 1.0 on successive refined meshes and for each mesh the corresponding 
error in L2 norm is computed as 



y {Qe{x,y,tf)-Whix,y,tf))^dxdy, (3.8) 



while the mesh size h{Cl{tf)) is taken to be the maximum diameter of the circumcircles 
of the triangles in the final domain Cl{tf). The resulting numerical convergence rates are 
listed in Table 1, whereas Figure 2 shows the temporal evolution of three different meshes 
used for the nimierical convergence study. 




Figure 2: Different mesh sizes used for the convergence rates study at different time outputs: t = (top row), 
t = l (middle row) and t = 2 (bottom row). The total number of elements Ne is increasing from the left grid 
(Ne=320), passing through the middle one (Ne = 1298), to the right one (Ne=5180). 
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3.2 Numerical Flux Comparison 

In order to study how the choice of the niunerical flux does affect the solution, we con- 
sider the well-known Sod shock tube problem. The rectangular shaped computational 
domain 0(0) = [— 0.5;0.5] x [— 0.1;0.1] contains a total number of elements of Ne = 18018. 
The initial condition reads 

0(x())-^ (1.0,0,0,1.0) if x<xo 

yix,uj-| (0.125,0,0,0.1) if x>xo ^"^-^^ 

with the position vector xq = (0,1/) which denotes the location of the discontinuity. 

The Sod problem involves physically a rarefaction wave traveling towards left and 
both a contact and a shock wave which are moving to the right side of the domain. Exact 
solution is known from the usage of an exact Riemann solver, while the numerical results 
plotted in Figure 3 have been collected with three different numerical fluxes, namely the 
Rusanov-type and the Osher-type flux given by (2.47) and (2.44), respectively. Further- 
more we use also an HLLC-type numerical flux, firstly introduced for moving meshes by 
Van der Vegt et al. [77, 78] in the DG finite element framework. For a detailed description 
of the HLLC flux and its extension to dynamic grid motion we refer to [77]. 

Figure 3 shows a comparison between the exact solution and the numerical results at 
time t = 0.25 obtained with a third order accurate scheme. The Rusanov flux is more dif- 
fusive if compared with the Osher and the HLLC fluxes, especially looking at the contact 
wave located at 0.23, which tends to be smoothed, while on the contrary it is almost 
perfectly resolved by the Osher- and the HLLC-ALE scheme. 

Convergence rate studies for the three different numerical fluxes are reported in Table 
2. As test problem we use again the isentropic vortex fully described in section 3.1, but 
here we consider only a numerical scheme with an order of accuracy of 03. The lowest 
L2—norm error is achieved by the very accurate Osher-type flux, as well as for the HLLC 
flux, whose error does not differ so much from the Osher scheme. The Rusanov type 
gives the highest error, as can be also observed in Figure 3. 

Table 2: Comparison of numerical convergence results for the compressible Euler equations using the third order 
version of the two-dimensional Lagrangian one-step WENO finite volume schemes and three different types of 
numerical fluxes (Rusanov, Osher and HLLC). The error norms refer to the variable p (density) at time t = 1.0. 





Rusanov 






Osher 






HLLC 




hmtf)) 




0(L2) 


h{Ci,tf) 




0{L2) 


h{Ci,if) 






3.61E-01 


1.076E-01 




3.28E-01 


1.614E-02 




3.31E-01 


1.818E-02 




2.51E-01 


2.315E-02 


4.2 


2.51E-01 


6.943E-03 


3.2 


2.51E-01 


7.897E-03 


3.0 


1.68E-01 


8.658E-03 


2.4 


1.68E-01 


2.290E-03 


2.7 


1.68E-01 


2.621E-03 


2.7 


1.28E-01 


3.950E-03 


2.9 


1.28E-01 


9.274E-04 


3.3 


1.28E-01 


1.068E-03 


3.3 



Looking at the computational time in Table 3, the Rusanov flux is the most efficent 
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scheme, but also the less accurate, while the Osher-type version gives the lowest error 
and is not the most expensive one, since the HLLC flux requires normally major com- 
putational efforts. We take Rusanov computational time as reference, hence rj evaluates 
simply the ratio between the current flux CPU time and the Rusanov one. Using HLLC 
one gets up to 7/ = 1.7, while at most 7/ = 1.5 is achieved with the Osher-type flux, which 
ensures the scheme to be the most accurate. 

Table 3: Computational time for the convergence studies and the Sod shock tube results obtained with three 
different numerical fluxes (Rusanov, Osher and HLLC). CPU Time is measured in seconds [s] and )/ denotes 
the ratio w.r.t. the Rusanov computational time. 





Rusanov 


Osher 




HLLC 






CPU time 


CPU time 


n 


CPU time 


V 




6.21E+00 


9.11E+00 


1.5 


1.06E+01 


1.7 


Isentropic 


2.22E+01 


2.49E+01 


1.1 


2.77E+01 


1.2 


vortex 


9.64E+01 


1.02E+02 


1.1 


1.12E+02 


1.1 




1.55E+02 


1.69E+02 


1.1 


1.93E+02 


1.2 


Sod problem 


1.22E+04 


1.44E+04 


1.2 


1.59E+04 


1.3 




Figure 3: Sod shock tube test. Third order accurate numerical results obtained using three different numerical 
fluxes (Rusanov, Osher and HLLC) and comparison with the exact solution (solid line) at time t = 0.25. 
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3.3 A Two-Dimensional Explosion Problem 

This test problem is in some sense a two-dimensional extension of the classical one- 
dimensional shock tube problems. The initial domain 0(0) is a circle of radius Ro=l and 
the initial condition is given by 

with = j-^+y^. The circle of radius R = 0.5 separates the inner state Qi from the outer 
state Qg. We assume 7 = 1.4 and use the initial condition reported in Table 4. 



Table 4: Initial condition for the explosion test 





P 


u 


V 


P 


Inner state (Q;) 


1.0 


0.0 


0.0 


1.0 


Outer state (Qo) 


0.125 


0.0 


0.0 


0.1 



In order to compute a reliable reference solution we assume rotational symmetry. 
Hence, one can simplify the multidimensional Euler equations to a one-dimensional sys- 
tem with geometric source terms, as proposed ia [75], which reads 

Qt+F(Q). = S(Q), (3.11) 

where 

Q=l puY F=| pS + p |, S = --( pu^ |. (3.12) 

\pE J \ uipE + p) J ' \ u{pE + p) J 

Here, r denotes the radial direction and u is the radial velocity, while a is a parameter 
that allows the system (3.11) to be equivalent to the one-, two- or three-dimensional 
Euler equations with rotational symmetry, according to its value: 

• a = 0: plain ID flow, 

• a = l: cylindrical symmetry (2D flow), 

• a = 2: spherical symmetry (3D flow). 

We choose oc — l for the 2D case and the inhomogeneous system of equations (3.11) is 
solved using a second order MUSCL scheme with the Rusanov flux on a one-dimensional 
mesh of 15000 points in the radial interval re [0;1]. This solution is assimied to be our 
reference solution, which will be used to verify the accuracy of the two-dimensional 
Lagrangian finite volume scheme. Figure 4 shows the comparison between the refer- 
ence solution and the numerical solution obtained with the third order version of our 
Lagrangian one-step WENO scheme for density and velocity at time t = 0.20: a circu- 
lar shock wave is travelling away from the center together with a contact wave, while a 
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circular rarefaction wave is running towards the origin. The contact wave is very well 
resolved due to the use of the little diffusive Osher-type flux. The mesh used contained 
68,324 triangles of characteristic mesh spacing /z = 1/100. A 3D view of the numerical 
solution as well as a very coarse version of the mesh are depicted in Fig. 5. 




0.2 0.4 0.6 0.8 1 "0 0.2 0.4 0.6 0.8 1 



Figure 4: Profiles along the positive x-axis of density (left) and velocity (right) at time f = 0.2. 




Figure 5: Left: density distribution at t = 0.20 for the explosion test. Right: coarse distorted mesh at i = 0.20. 

3.4 The Kidder Problem 

The Kidder problem consists in an isentropic compression of a shell filled with an ideal 
gas. For this problem, an exact analytical solution has been proposed by Kidder in [50]. 



21 



This is a very useful and widely used test [10, 56], to assure that a Lagrangian scheme 
does not produce spurious entropy. The shell has a time-dependent internal radius ri{t) 
and an external radius re{t), while r denotes the general radial coordinate with ri(t) < 
r < re{t). The initial values for the internal and external radius are r,(0) =r, o = 0.9 and 
re(0) = fg^o = 1-0, respectively. In this test problem the ratio of specific heats is 7 = 2 and 
the initial density distribution is given by 



po = p{r,0) = 




(3.13) 



with Pi,o =1 and pe,o =2, which are the initial values of density defined at the internal and 



at the external frontier of the shell, respectively. The initial entropy sq = 
so that the initial pressure distribution is given by 

poir)=sopoiry. 



Po 

'Pl 



- 1 is uniform. 



(3.14) 



Initially the fluid is at rest, hence u = v = Q. 

According to [50] we look for a solution of the the form R{r,t) = h{t)r, where R{r,t) 
denotes the radius at time t > of a fluid particle initially located at radius r. Therefore 
the self-similar analytical solution for t G [0,t] reads 



p{R{r,t),t) = h{ty 

Ur{R{r,t),t) = j^hit) 

piR{r,t),t) = h{t) 
where the homothety rate is 

h{t)- 

with T representing the f ocalisation time 



Po 



h{t) 



R{r,t) 
h{t) 



27 
7-1 



Po 



R{r,t) 
hit) 



T = 



\ 



7-1 (.'lo-i,o) 



'-6,0 ^i,0 



and the internal and external soimd speeds c, and Cg are defined as 



7- 



Pi 



Pe' 



7 



(3.15) 



(3.16) 



(3.17) 



(3.18) 



The boundary conditions are imposed at the internal and the external frontier of the 
shell, where a space-time dependent state is prescribed according to the exact analytical 
solution (3.15). 
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As done by Carre et al. in [10] the final time of the simulation is chosen to be = ^ t, 
so that the resulting compression rate is h{tf) =0.5. In this way the exact solution is given 
by a shell bounded with 0.45 <R< 0.5. 




Figure 6: Left: Numerical solution at J = 0.4 and at t = tf. Right: Evolution of the internal and external radius 
of the shell and comparison between analytical and numerical solution. 

We use the Osher-type numerical flux (2.47) to obtain the results shown in Figure 6, 
where we can observe an excellent agreement between the analytical and the numerical 
solution. The absolute error \err\ concerning the location of the internal and external 
radius at the end of the simulation has also be computed and is reported in Table 5. 





^ ex 


^num 


\err\ 


Internal radius 


0.4500000 


0.4499936 


6.40E-06 


External radius 


0.5000000 


0.4999922 


7.80E-06 



Table 5: Absolute error for the internal and external radius location between exact R^x and numerical R„ 
solution. 



3.5 The Saltzman Problem 

A strong one-dimensional shock wave driven by a piston constitutes this test problem. 
As proposed in [53,56], the initial two-dimensional domain is a box 0(0) = [0;1] x [0;0.1], 
filled with an ideal gas. The piston is pushing the gas from the left to the right and it 
is moving with constant velocity Vp = (1,0). The computational domain is covered by 
a distorted unstructured mesh, composed of 100 x 10 elements along the boundaries, as 
reported in Figure 7. 
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Time t = 0.0 




Figure 7: The initial and the final mesh for the Saltzman test problem. 

The initial condition is given by an ideal gas at rest (vq = 0) with a ratio of specific 
heats 7 = |, an initial density of po = and an internal energy of Cq = 10^^, corresponding 
to an initial pressure of po = 6.67 -IQ"^ . 

As boundary conditions we set fixed slip wall boundaries on the upper and lower 
boundary outflow on the right and a moving slip wall boundary condition for the piston 
on the left. According to [53] we use initially a Courant number of CFL = 0.01, in order 
to respect the geometric CFL condition and to avoid mesh elements crossing over. After 
this initial phase, the CFL condition is raised to its usual value of CFL = 0.5 in two space 
dimensions. 

The exact solution is a one-dimensional infinite strength shock wave and it can be 
computed by solving the Riemann problem given in Table 6. 





Left state 


Right state 


p 


1.0 


1.0 


u 


1.0 


-1.0 


V 


0.0 


0.0 


V 


6.67-10-7 


6.67-10-7 



Table 6: One-dimensional Riemann problem for obtaining the exact solution of the Saltzman problem. 

The details of the algorithm that computes the exact solution of the Riemann problem 
are given in the book of Toro [75]. The exact solution has then to be shifted by a certain 
quantity d to the right, corresponding to the movement of the piston during the time of 
the simulation tf, i.e. 

d = Up-tf. (3.19) 

The exact post shock density is pe=4.0 and the final time is chosen to be tf =0.6, according 
to [53]. Therefore, the exact final shock location is at x = 0.3. 

We use a Rusanov-type flux (2.44) and the solution is computed, both, for a third 
order and a fourth order scheme. We can notice a very good agreement regarding the 
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final density distribution, as depicted in Figure 8, and the solution for the velocity, where 
we observe a very sharp discontinuity at the shock location, as clearly shown in Figure 
9. The decrease of the density near the piston in Figure 8 is due to the well known wall- 
heating problem, see [74]. One can notice from Figure 10 the progressive mesh compression 
at the inflow, where the piston is pushing, and the shock wave that is travelling through 
the domain, reaching its final location exactly according to the reference solution {x=0.3). 



0.2 0.3 
X 



£2.5 

2 



Figure 8: Numerical results for the density distribution of the Saltzman problem: third order scheme (left) and 
fourth order scheme (right). Both results are compared with the analytical solution, obtained by solving the 
Riemann problem described in Table 6. 




Figure 9: 3D visualization for the velocity distribution obtained by a third order finite volume scheme (left) and 
comparison with the analytical solution (right). 
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Figure 10: Progressive evolution of the numerical solution of the Saltzman problem (density distribution): the 
mesh is compressed by the piston on the left, while a shock wave is travelling towards the right. The final 
location coincides perfectly well with the exact solution. 
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3.6 Single Mach Reflection Problem 

In this section we solve the single Mach reflection problem. It consists of a shock wave 
traveling at a shock Mach number of Ms =1.7 that hits a wedge of an angle tt=25°. Experi- 
mental reference data in form of Schlieren images as well as numerical reference solutions 
are shown, e.g. in [75]. The upstream density and pressure are po=l and pQ = l/j, respec- 
tively, where the ratio of specific heats is 7 = 1.4. In the following, the indices o and i will 
denote upstream and downstream states, respectively. Using the Rankine-Hugoniot con- 
ditions, we setup the computation with a shock wave initially centered at x = —0.04 that 
travels at Mg = 1.7 to the right into a medium at rest. The wedge is defined by the triangle 
W25, composed of the vertices (0,0), (3.0,0) and (3.0,1.398923). The initial computational 
domain Q(0) = [-2;3] x [0;2]\yV25 is discretized with 177,440 triangles of characteristic 
size h = 1/100 and a P0P2 WENO finite volume scheme is used. At the final time t — 1.2 
the exact shock location must he x = 2. The upper and lower boundaries of the computa- 
tional domain are discretized by slip walls. At the other boundaries Dirichlet boundary 
conditions consistent with the initial condition are imposed. 

The results of the computation are depicted in Fig. 11, where the density contour lines are 
shown, as well as a zoom into the final mesh, which is very distorted along the slip line. 
This is to verify that our proposed high order unstructured finite volimie algorithm can 
handle reasonably strong mesh deformations, which are a common feature of Lagrangian 
computations. For even stronger mesh deformations, as they occur in the double Mach 
reflection problem, a proper remeshing and projection strategy will be implemented in 
the future. The general flow features depicted in Fig. 11 agree very well with the compu- 
tations and the experimental data shown in [75]. 

3.7 Two-dimensional Riemann Problems 

In this section we solve a set of two-dimensional Riemann problems, whose initial con- 
ditions are given by 



A large set of such two-dimensional Riemann problems has been presented in great de- 
tail in the paper by Kurganov and Tadmor [51]. The initial conditions for the three config- 
urations presented in this article are listed in Table 7. The initial computational domain 
is defined as 0(0) = [— 0.5;0.5] x [— 0.5;0.5]. The Lagrangian simulations are carried out 
with a third order one-step WENO scheme using an imstructured triangular mesh com- 
posed of 90,080 elements with an initial characteristic mesh spacing of h = 1/200. The 
reference solution is computed with a high order Eulerian one-step WENO finite vol- 
ume scheme as presented in [25, 27, 28], using a very fine mesh composed of 2,277,668 
triangles with characteristic mesh spacing h = 1/1000. In all cases, the Rusanov flux has 




(3.20) 
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been used. The exact solutions of the one-dimensional Riemann problems are imposed 
as boundary conditions on the four boundaries of the domain. The results obtained with 
the Lagrangian-type scheme together with the final mesh and the Eulerian reference so- 
lution are depicted in Figures 12 - 13. We can note a very good qualitative agreement of 
the Lagrangian solution with the reference solution, as well as with the results published 
in [51]. 

Table 7: Initial conditions for the two-dimensional Riemann problems. 



Problem RPl (Quadruple Sod problem) 





x<0 


x>0 




p U V p 


p u V p 


y>0 
y<0 


1.0 0.0 0.0 1.0 
0.125 0.0 0.0 0.1 


0.125 0.0 0.0 0.1 
1.0 0.0 0.0 1.0 


Problem RP2 (Configuration 2 in [51]) 




x<0 


x>0 




p 11 V p 


p 11 V p 


y>0 
y<0 


0.5197 -0.7259 0.0 0.4 
1.0 -0.7259 -0.7259 1.0 


1.0 0.0 0.0 1.0 
0.5197 0.0 -0.7259 0.4 


Problem RP3 (Configuration 4 in [51]) 




x<0 


x>Q 




p U V p 


p U V p 


y>0 
y<0 


0.5065 0.8939 0.0 0.35 
1.1 0.8939 0.8939 1.1 


1.1 0.0 0.0 1.1 
0.5065 0.0 0.8939 0.35 



4 Conclusions 

We have developed a new high order two-dimensional Arbitrary-Lagrangian-Eulerian 
one-step WENO finite volume scheme on unstructured triangular meshes. The algo- 
rithm works for general hyperbolic balance laws with non stiff algebraic source terms. 
Several smooth and non-smooth test problems with shock waves, contact waves, shear 
waves and rarefactions have been simulated. The results have been compared with exact 
or numerical reference solutions in order to validate our approach. The algorithm was 
found to work properly in all cases and to be robust and accurate, too. The accuracy has 
been verified empirically via numerical convergence studies on a smooth test case with 
exact solution. Further work will concern the improvement of the presented algorithm 
in order to deal with stiff algebraic source terms, which is straightforward, see the work 
presented in [30]. Moreover, we plan to extend our scheme to three space dimensions in 
the more general framework of the new P^Pm method proposed in [25], where one can 
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Figure 12: Top left: 50 equidistant density contour lines from |O = 0.14 to ,0 = 0.98 for the numerical solution 
obtained with our third order unstructured Lagrangian one-step WENO finite volume scheme for the two- 
dimensional Riemann problem RPl at t = 0.2. Top right: Same contour lines for the reference solution. Bottom: 
computational mesh of the Lagrangian scheme at the final time t = 0.2. 
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0.5 
0.4 
0.3 
0.2 
0.1 


-0.1 ^ 
-0.2 
-0.3 
-0.4 
-0.5 
-0.6 
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X 

Figure 13: Top left: 51 equidistant density contour lines from p = 0.27 to ,0 = 0.98 for the numerical solution 
obtained with our third order unstructured Lagrangian one-step WENO finite volume scheme for the two- 
dimensional Riemann problem RP2 at t = 0.2. Top right: Same contour lines for the reference solution. Bottom: 
computational mesh of the Lagrangian scheme at the final time f = 0.2. 
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Figure 14: Top left: 51 equidistant density contour lines from |O = 0.53 to p = 1.9 for the numerical solution ob- 
tained with a third order unstructured Lagrangian one-step WENO finite volume scheme for the two-dimensional 
Riemann problem RP3 at f = 0.245. Top right: Same contour lines for the reference solution. Bottom: compu- 
tational mesh of the Lagrangian scheme at the final time t = 0.2i5. 



32 



deal with pure finite volume or pure finite element methods, or with a hybridization of 
both. Further work will also consist in a generalization to moving curved meshes as well 
as to hyperbolic PDE with non-conservative products. 
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